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ABSTRACT 


It  is  shown  how  to  construct  a modified  version  SUB!^  of  a (presumably 


efficient)  subroutine  SUB^  for  solving  the  linear  system  A^x  = b , i=l, 


fk, 


so  that  the  linear  system 

(A^  ® ® A^)x  = b 

can  be  solved  by  just  one  call  to  each  of  the  routines  SUBj^  , i=l, ,k.  Poly- 

nomial interpolation  and  spline  interpolation  in  several  variables  are  given  as 
examples . 
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SIGNIFICANCE  AND  EXPLANATION 

The  tensor  product  of  linear  approximation  schemes  is  a very  convenient 
matiiematical  construct  to  make  efficient  multivariate  approximation  schemes  out 
of  univariate  ones.  For  example,  we  might  know  that  the  coefficient  vector  a 
for  the  polynomial 

p(x)  = y a.  p. (x) 

‘"11 

1 

which  agrees  with  the  function  f at  t,,...,t  can  be  computed  as 

In 

a = (f (t. ) ) . 


Then  we  know  that 


= I I a^.P^(x)p  (y) 


1 : 

is  a polynomial  which  agrees  with  the  function  f at  all  the  points  (t.,u.)  of 

a rectangular  grid  provided  the  coefficient  matrix  is  computed  from  the 

data  matrix  (f(t.,u.))  by 
1 1 

U.  .)  = ) (f(t.,u.))  :=  (f(t..u.))(B  )'^  . 

1]  11  ^.ilE 

Such  a construction  can  also  be  used  to  obtain  interpolants  to  functions  of  three 
or  more  variables,  but  explicit  expressions  for 

(B^®B  ®...®  B ) (f  (t.  ,u.  , . . . ,z  )) 

i E £ 1 1 

in  terms  of  ordinary  matrix  multiplication  are  harder  to  come  by  and  implement. 

The  report  proposes  a simple  way  of  modifying  a subroutine  for  calculating  the 
coefficient  vector  a from  univariate  data  so  that  the  construction  of  a tensor 
product  interpolant  reduces  to  calling  these  modified  univariate  routines  in 
sequence.  Somewhat  surprisingly,  this  idea  seems  not  to  have  been  noticed  before. 

It  makes  the  construction  and  evaluation  of  tensor  product  approximants  an  easy  thing. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive  summary 
lies  with  MRC,  and  not  with  the  author  of  this  report. 


EFFICIENT  computer  MANIPULATION 


OF  TENSOR  PRODUCTS 

Carl  de  Boor 

In  [3],  V.  Poreyra  and  G.  Scherer  discuss  the  numerical  solution  of  a linear  systejn  of 
the  form 

(A,  S'  . . . S'  A,  )x  = b ( 1 ) 

1 k = = 

with  A^  an  invertible  matrix  of  order  n^,  i=l,...,k,  and,  correspondingly,  both  x and 
b k-dimensional  arrays,  of  size  n^^  x x...xnj^.  Such  systems  arise  naturally  when  form- 
ing tensor  products  of  univariate  interpolation  schemes. 

Pereyra  and  Scherer  propose  to  store  arrays  such  as  x and  b with  the  last  index 
running  fastest  and  then  have  a scheme  of  applying  down  to  and  includ- 

ing 1^2^'  appropriately  restoring  the  intermediate  information  so  that  application  of 
A^^  involves  only  repeated  ordinary  matrix  multiplication  to  a vector  stored  in  consecutive 
locations  in  memory.  When,  as  is  more  reasonaible,  application  of  rather  than  of 

A^^  is  wanted,  with  L^U^  a triangular  factorization  for  A^,  a further  complication 
arises  and  is  dealt  with. 

It  is  the  purpose  of  this  note  to  describe  a different  procedure  which  I have  used  for 
some  time  and  which  is  more  direct  and  simpler  than  the  Pereyra-Scherer  procedure  appears  to 
be. 

We  assume  that,  for  each  i,  we  have  available  a Fortran  subroutine 

SUB. (b,n,x) 

1 

which  solves  the  i-th  linear  system  A^x  = b (of  order  n - n^)  for  x,  given  b.  Pre- 
sumably, the  routine  does  this  in  an  efficient  way,  taking  advantage  of  any  special 
structure  A^  might  have  such  as  bandedness,  positive  definiteness  etc. 
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We  further  assume  that  the  k-dimensional  arrays  x and 
Fortran  fashion,  i.e., 


b are  (to  be)  stored  in 

+n.  ,(i  -1)..))) 

lc-1  k 


if  we  refer  to  x also  as  an  equivalent  one-dimensional  array. 

The  following  simple  procedure  will  then  lead  to  an  efficient  way  for  solving  (1)  : For 
each  i,  enlarge  the  subroutine  SUB^  to  a subroutine 


sub: (b,n,m,x) 
1 


which  solves  simultaneously  A^x  = b for  m given  right  sides  b(",l),  b ( • , 2) , . . . ,b( • ,ra) , 
each  of  length  n = n^,  and  stores  the  corresponding  solutions  in  x(l,*),  x(2, .x(m, *) . 
Thus,  the  dimension  statement  for  the  arguments  b and  x in  SUB^  reads 

DIMENSION  b(n,m),  x(m,n) 

and  the  change  otherwise  consists  in  putting  every  statement  involving  b or  x appropriate- 
ly into  a DO  loop.  In  this,  care  should  be  talten  to  leave  statements  which  do  not  depend  on 
the  particular  right  side  outside  such  loops. 

Many  library  routines  for  specific  linear  problems  already  provide  this  facility  for 
dealing  with  several  right  sides  in  one  call,  since  t)ie  work  in  solving  A^x  = b for  an 
additional  right  side  b is  usually  much  less  than  the  work  for  solving  such  a system  the 
first  time.  But  such  routines  return  the  solution  corresponding  to  the  j-th  column  b(*,j) 
of  the  input  array  customarily  in  the  j-th  column  of  the  output  array  x and  not,  as  I pro- 
pose here,  in  the  j-th  row. 

Lemma . For  i=l,...,k,  let  SUB^  be  an  expanded  version,  as  described,  of  the  routine 
SUB^  for  solving  A^^x  = b,  and  set 


N :=  n. 


Then,  the  following  statements 


CALL  SUB|(b^^,  n^,  N/n^ , b^) 


CALL  SUB2(bj,  n^.  N/n^/  b^) 


CALL  SUB'(b  n , N/n  , b ) 

^ K X iC 


S ==  6k 


will  produce  the  solution  x of  (1)  . 


Proof.  Let  be  the  k-dimensional  array 


X.  :=  {a"^  ® . . . ® aT^  ® 1 ® . . . ® l)b 

= 1 1 1 = 


i=0, . . . ,k.  Then 


= i‘^l ^i-1'  ■ '^i+1 ^k’ 


-1 


(2) 


= ^ ^i-r  • '^i+1 3k> 


and  our  assertion  is  proved  if  we  can  establish  that  b^^  = We  prove  this  by  showing 


that,  for  all  i,  b^  as  generated  by  the  succession  of  calls  above  is  related  to  x.  in 


the  following  way: 

If  is  interpreted  as  a k-dimensional  Fortran  array,  of  dimension 


*"i+l"  ■ ■ '"k'"!  ' ■ ■ ■ '"i' ' 


Si<^i+i \'^i 3i>  = 5i<^i V'  ^ i ' 


(3) 


for  i=0,...,k.  For  i = k,  (3)  is  indeed  the  desired  statement  that  b = x 
=k  =k 


Now,  (3)  holds  for  i=0  because  of  the  initial  assignment  b^  :=  b.  Assuming  (3)  to 


hold  for  i<v,  we  consider  the  action  of  the 


CALL  SUB’ (b  , ,n  ,N/n  ,b  ) . 

V =V-1  V V =v 


SUB^  considers  b^  to  be  a two-dimensional  array,  b say,  of  dimension  (n  ,N/n^)  . 


Thus,  with 


we  have 


:=  1 , ,-i  +...+n  (j,-l 

v+1  v+l  v+2  k 1 


b(.,s)  = j^.,) 


(4) 


by  induction  hypothesis.  SUB^  then  applies  to  each  of  these  m - N/n^  n^-vectors 

b(-,s),  thus  obtaining  the  corresponding  n^-vector 


X ( j,  , 
=v  1 


^v-1' 


-■v+l 


by  (2)  and  (4) , and  stores  this  vector  in 


= b (s  + 

(N/n  ) ( • -D) 

~v 

V 

= fev<Vl 

. . ) + (N/n^) ( ■ 

■ -1) 

= fev'^v+l 

^k'^i Vi'  • > 

which  proves  (3)  for  i = v and  so  advances  the  induction  hypothesis;  Q.E.D. 

We  introduced  the  auxiliary  arrays  only  for  argument's  sake.  In  calculations,  two  arravs, 
say  and  b^,  are  sufficient,  with  bj  serving  in  place  of  all  b^  with  i odd,  and 

b^  serving  for  all  the  others. 


Also,  in  typical  situations,  the  various  subroutines  SUB^  are,  in  fact, 

just  one  routine  called  with  additional  arguments  which  differ  with  i.  In  such  a case,  only 
one  extended  version  has  to  be  written. 

Finally,  we  put  the  above  discussion  in  terms  of  solving  a linear  system,  i.e.,  in  terms 
of  premultiplying  a given  vector  by  the  inverse  of  a given  matrix.  We  did  this  in  order  to 
make  the  point  that  we  do  not  require  the  matrix  by  which  we  wish  to  premultiply  to  be  pre- 
sent explicitly.  Any  Fortran  subprogram  SUB^(b,x)  which  has  the  effect  of  forming 
X = Bj^b  for  given  b can  serve  as  a basis  for  an  extended  version  SUBMb,n,ra,x)  suitable 
for  the  calculation  of  (B^  ® matrices  need  not  be  square.  We  state 

this  slight  extension  of  the  Lemma  as  a corollary  for  the  record. 
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ol  lar\'.  For  i*l 


k,  let  B.  be  a (n./r  )-inatriX/  and  let  SUB  ! {b , n . ,m,  x , r ) 
i 1 j 1 1 1 

bo  a subroutine  which,  for  j=l m,  forms  the  r ^"Vector  b^b(*,j)  (in  some  mariri<-rj 

from  the  n ^-vector  b{*,j),  and  stores  it  in  x(j,*).  Then,  th€»  following  statements 


CALL  SUB|(bp,  n^,  m,  , r^) 

m :=  m*r^/n2 

CALL  SUB^(b^,  m,  b^,  r^) 

m :=  m*r2/n^ 

CALL  sub;  (b  , , n,  , m,  b,  , r,  ) 
k =.<-1  k =k  k 


form  the  k-dimensional  array  x=  (B^®  ...®Bj^)b. 

It  is  not  even  necessary  that  be  a matrix,  i.e.,  a two-dimensional  array.  The  more 

general  situation  in  which  B^  is  a linear  map  which  associates  s ^-dimensional  arrays  with 
t^-dimensional  arrays  is  covered  by  the  corollary  as  well  since  we  can  always  interpret  such 
s^-dimensional  and  t^-dimensional  arrays  Fortran  fashion  as  equivalent  one-dimensional  arrays. 

We  give  some  simple  examples  in  the  next  section. 

Tensor _pTO^cts  of  univariate  interpolation  schemes.  The  following  material  concerning 

tensor  products  of  univariate  interpolation  schemes  is  well  known  and  is  mentioned  here  only 
in  order  to  illustrate  the  use  and  usefulness  of  the  simple  idea  expounded  earlier.  (A 
simple  account  giving  proofs  and  details  can  be  found,  e.g.,  in  [1].) 

The  construction  of  a (univariate)  linear  interpolant  g to  some  function  f usually 
involves  the  calculation  of  the  coefficients  a = (a.)  in  a representation 


B 


g = y a.sP  . 

r 1 i 


1 

for  the  interpolant  from  certain  information  (A^f)  about  f.  Here,  each  A^  is  a linear 
functional , e.g., 

(r . ) 

A.f  = f(x.)  or  A-f  = f (x.)  or  A-f  = [ ij,.  (x)  f (x) dx  etc. 

1 1 1 1 1 J 

and  g is  so  constructed  that 

A^g  = A^f  , all  i 

At  the  level  of  the  present  discussion,  there  is  no  reason  to  require  the  representation  for 
g to  be  irredundant,  i.e.,  to  require  the  sequence  ('A^)  to  be  linearly  independent.  All 
that  is  necessary  is  the  assumption  that 

a = B{A.f) 

1 

for  some  matrix  B.  The  matrix  B is  commonly  not  known  explicitly  (although  it  could,  of 
course,  be  determined).  Rather,  some  procedure  or  subprogram  SUB  is  available  which  trans- 
forms the  vector  (A^f)  of  data  appropriately  into  the  vector  a of  coefficients. 

For  example,  consider  the  construction  of  the  polynomial  P = p^  of  degree  < n which 

agrees  with  f at  the  n distinct  points  x, , . . . ,x  . In  its  Newton  form,  p,  looks  like 

In  r 


II  11 

P,(x)  = I lx  , ...,x  If  • 1 I (x  - X ) 
f i=l  ^ " j=i+l  ^ 


(5) 


with  the  coefficient  Ix.,...,x  ]f  the  so-called  divided  difference  for  f at  the  points 

i n 

x.,...,x  , i=l,...,n,  i.e., 

1 n 


f (x. ) 
1 


i = 1 


[x ,x.]  f :=  ( (6) 

" ^ I 

( lx  , . . . ,x.  1 f - lx.  , — ,x.  1 f)/(x  .-X. ) , i < j 

Is.  1+1  1 1 1-1  11 


These  coefficients  can  therefore  be  determined  as  final  entries  in  a so-called  divided 
difference  table,  for  instance  as  in  the  following  subprogram 
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SlJBKiHJTINE  PnUNT  (X,  F,  M) 
DIMENSION  X(N) ,F(N) 


NMl  = N-1 

IF  (NMl  .LE.  0)  RETURN 

DO  10  K=1,NM1 
NMK  = N-K 
DO  10  1=2, NMK 

10  F(I)  = (F(I+1)  - F(I) )/(X(I+K)  - X(I)) 

END  RETURN 

Here,  the  array  F contains  F(i)  = f(x^),  i=l,...,n,  on  input  and  F(i)  = , • • • >x^l f , 

i=l,...,n,  on  output.  (For  details  concerning  divided  differences  and  the  Newton  form  (5), 
see,  e.q. , 12] . ) 

Once  the  coefficient  vector  a in  the  representation  interpolant  q 

has  been  determined,  one  may  evaluate  g in  various  ways.  Typically,  one  then  wants  to  find 
Xg  for  various  linear  functionals  X such  as  Xg  = g(x),  some  x,  or  Xg  = g*^^(x),  or 
Xg  = / iig  for  some  '{’>  etc.  . All  of  these  values  can  be  obtained  from  the  vector 
a = (a^)  by  applying  to  it  a matrix  consisting  of  just  one  row,  viz.  the  matrix 
(Xip^,  ^'^2'  Thus  evaluation  of  the  interpolant  at  some  linear  functional  X is  just 

another  linear  procedure  or  subprociram  which  applies  some  matrix  B to  the  vector  a. 

For  example,  the  evaluation  of  the  interpolating  polynomial  (5)  at  some  [xjint  x = ARC 
proceeds  customarily  by  Nested  Multiplication,  as  in  the  following  function  subprogram 


10 


FUNCTION  POLVAL  (X,  F,  N,  ARC) 

DIMENSION  X(N),  F(N) 

POLVAL  = F(l) 

IF  (N  .LE.  1)  RETURN 

DO  10  K=1,N 

POLVAL  = POLVAL* ( ARC- X(K))  + F(K) 

RETURN 

END 
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Note  that,  once  aqain,  the  matrix  B to  be  applied  to  the  coefficient  vector  a fm  the 
array  F)  is  not  formed  exi>licitlv. 


Suppose  now  that  we  have,  for  each  of  the  k independent  variables  tj,...,t^,  a 


linear  interpolation  scheme.  This  means  that,  for  r=l,,..,k,  we  have  a matrix  B which 

r 

associates  with  each  data  vector  (X^f)  a coefficient  vector  (a^)  = B { ^Tf) , qivinq  the 

1 1 r i ^ 

interpolant  ^ for  f = f(t^)  , Further,  for  all  appropriate  integer  vectors 

_i  = ( i^ , . . . , ij^)  , let  be  a linear  functional  on  some  appropriate  class  of  functions  f 


of  k variables  for  which 


\.f  = (A^  f,)(A^  f^) 

1 1,  1 1-  2 i-,  ^ 

= 12  k 


whenever 


h \ • 

b 

For  example,  if  k = 3 and  xff  = f(a  ),  X^'f  = f"(B  ) and  X^f  = f ^f(t)dt,  then 

1 r 2 r 3 a 


(1,2,3)^  J^a,  ^*“l'®2''^3'‘^’^3 


^(2,2,1)^  ''  f(6j^. 62.03) 


would  serve.  Also,  let 


==  ^,,l'^l’^,,2<^2’•••'^i  ,k<^’  • 

= 12k 


Then  we  can  construct  an  interpolant 


g - y a.^  . 

h 11 


for  a function  f of  the  k variables  t,,...,t.  as  follows:  Calculate  the  k-dimens ional 

1 k 


array  a = (a^)  as 


a = (B,  ® B.  ) (X.f) 

= 1 k 1 


from  the  k-diitiensional  array  (\^f)  of  data.  This  function  g is  then  indeed  an  inter- 
polant  to  f in  the  sense  that 

i,g=i.f,  all  i 
1 = 

The  calculation  of  the  coefficient  array  a is,  of  course,  easily  effected  as  described  in 
the  corollary  above. 

To  follow  up  on  the  example  of  polynomial  interpolation,  an  appropriately  extended  ver- 
sion POLNTE  of  the  subprogram  POLINT  would  require  a separate  output  array,  D say,  for 
the  calculated  divided  differences.  Otherwise,  onlv  the  statement  labelled  10, 

10  F(I)  = (F(I+1)  - F(I) )/(X(IPK)  - X(I)) 

needs  to  be  put  into  an  additional  loop  over  the  data  sets,  with  the  difference  X(IPK)  - X(I) 
calculated  outside  that  loop,  of  course.  We  get 
SUBROUTINE  POLNTE  (X,  F,  N,  M,  D) 

DIMENSION  X(N) ,F(N,M) ,D(M,N) 

DO  5 1=1, N 
DO  5 J=1  ,M 

S D(J,I)  = F(I,J) 

NMl  = N-1 

IF  (NMl  .LE.  0)  RETURN 

DO  10  K=1,NM1 
NMK  = N-K 
DO  10  1=1, NMK 

DIFF  = X(I+K)  - X(I) 

DO  10  J=1,M 

10  D(J,I)  = (D(J,I+1)-D(J,I) )/DIFF 

RETURN 

END 


-9- 


Note  that  this  routine  functions  appropriately  even  for  M = 1,  the  only  difference  com- 
pared to  POLINT  being  that  the  output  is  now  to  be  found  in  D and  not  in  F.  Mote 

further  that  it  takes  N(N-l)/2  adds  and  divides  per  data  set  to  form  3(X.f).  Since  the 

1 

matrix  B ^ is  upper  triangular  in  this  case,  explicit  application  of  B by  backsubstitu- 
tion  would  take  no  fewer  operations  and  would  require  the  generation  and  storage  of  B (or 
its  inverse) . 

Now,  to  illustrate  the  lemma  and  its  corollary,  suppose  that  we  require  the  polynomial 
interpolant  p = p(x,y,z)  to  data 

f(x.,y.,z,)  , i=l,...,h  j j=l,...,n  ; k=l,...,n 
1 D k X y z 

We  load  f(x.,y.,z  ) into  F(i,j,k),  x.  into  X(i) , y.  into  y(j)  and  z into  Z(k) , 
13^  ^ 3 k 

for  all  appropriate  i,j,k.  Then 

N :=  n *n  *n 
X y z 

CALL  POLNTE  (X,  F,  n , N/n  , D) 

X X 

CALL  POLNTE  (Y,  D,  n , N/n  , F) 

y y 

CALL  POLNTE  (Z,  F,  n , N/n  , D) 
z z 

to  get  the  appropriate  polynomial  coefficients  of  the  polynomial  interpolant  p into  the 
3-dimensional  array  D. 

If  we  wish  to  evaluate  this  interpolant  at  some  point  (x,y,z) , we  have  to  procure  an 
extended  version  of  the  function  routine  POLVAL-  The  output  for  such  a routine  will  consist 
now  of  more  than  one  number,  wc  must  therefore  give  up  on  having  a function.  Otherwise,  it 
is  again  only  the  assignment  statement  POLVAL  = F(l)  and  statement  10  which  need  to  be  put 
into  a loop  over  the  data  sets.  Here  is  an  extended  version  POLVLE  of  POLVAL. 
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SUBROUTINE  POLVLE  (X,  D,  N,  M,  ARC,  VALUE) 
DIMENSION  X(N) ,D(N,M) .VALUE(M) 

DO  5 J=1,M 

5 VALUE(J)  = D(1,J) 

IF  (N  .LE.  1)  RETURN 

DO  10  K=2,N 

FACTOR  = ARC  - X(K) 

DO  10  J=1,M 

10  VALUE  (J)  = VALUE  (J) ‘FACTOR  + D(X,J) 

RETURN 


END 

Now,  to  find  p(x,y,z)  , 

CALL  POLVLE  (X,  D,  n , N/n  , x,  TEMPI) 

X X 

CALL  POLVLE  (Y,  TEMPI,  n , n , y,  TEMP2) 

y z 

CALL  POLVLE  (Z,  TEMP 2 , n , 1,  z,  ANSWER) 

z 

to  get  p(x,y,z)  = ANSWER.  Note  that  TEMPI  must  be  of  size  n *n  and  contains  the 

y z 

necessary  information  to  evaluate  the  bivariate  polynomial  p(x,y,z)  for  any  choice  of  y 

and  z.  Again,  TEMP2  is  of  size  n and  contains  the  appropriate  coefficients  of  the 

z 

polynomial  p(x,y,z)  in  the  single  variable  z.  In  particular,  if  p is  to  evaluated  at 
all  points  of  a regular  grid,  it  is  most  efficient  to  evaluate  p along  lines  parallel  to 
the  z-cixis. 

As  an  example  of  some  of  the  difficulties  one  might  encounter,  we  now  discuss  briefly 
osculatory  polynomial  interpolation.  Here,  the  interpolant  is  again  of  the  form  (5),  but 
now  some  of  the  interpolation  points  Xj^,...,x^  might  coincide.  This  requires  an  exten- 
sion of  (6)  which  reads  as  follows: 


X.  = X.  implies  x.  =x.  ,=...=x.  , 

1 j ' 1 1+1  j 

(6)  and  (6a)  cover  all  eventualities.  The  point  of  this  extension  is  that  now  p^  aarees 
with  f in  the  sense  that 

p^'^'iz)  - f*'^*(Z)  in  case  the  number  z appears  (at  least)  r + ] 

times  in  the  sequence  x ,x 

1 n 

This  explains  the  term  "osculatory “ . 

The  following  program  for  the  construction  of  the  coefficients  in  (5)  is  based  on  (6) 
and  (6a)  and  can  be  found,  in  somewhat  different  notation,  in  12], 

SUBROUTINE  POLOSC  (X,  F,  N) 

C INPUT  MUST  SATISFY  THE  FOLLOWING. 

C IF  X(I-l)  .NE.  X(I)  = X(I+J)  .NE.  X(I+J+1),  THEN 
C X(I+L)  = X(I)  AND  F(I+L)  = (D**L) F( X ( I) ) , L=0,...,J  . 

C (HERE,  X(0)  , X(N+1)  .NE.  X(I),  1=1 N,  BY  DEFINITION.) 

DIMENSION  X(N) ,F(N) 

NMl  = N-1 

IF  (NMl)  .LE.  0)  RETURN 

DO  10  K=1,NM1 
FLOATK  = K 
NMK  = N-K 
FLAST  = F(l) 

DO  9 1=1, NMK 

DX  = X(I+K)  - X(I) 

IF  (DX  .EO.  0. ) GO  TO  7 

F(I)  = (F(I+1)  - FLAST)/DX 
FLAST  = F(I+1) 

GO  TO  9 

7 F(I)  = F(I+1) /FLOATK 

9 CONTINUE 

10  F(NMK+1)  = FLAST 

RETURN 

END 

The  construction  of  an  efficient  extension  of  POLOSC  is  made  difficult  by  the  fact 
that  the  local  variable  FLAST  depends  on  the  data  F but  is  active  through  various 
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statemt'iits  whicji  are  independent  of  the  data  F and  should  therefore  not  be  put  inside  a 
loop  over  the  various  data  sets.  One  wav  out  is  to  make  FLAST  an  array  of  length  .M, 
either  local  or  as  an  argument,  which  then  requires  the  four  groups  of  statements 
Fl-AST  = F(l) 

F(I)  = (F(I+1)  - F(I))/DX  ; FLAST  = F(I+1) 

F(I)  = F(I  + 1) /FUIATK 
F(NMK+1)  = FLAST 


each  be  put  into  a loop  over  the  different  data  sets. 

An  alternative  way  consists  in  a reorganization  of  the  entire  calculation  which  avoids 
the  temporary  saving  of  terms  which  depend  on  F,  possibly  at  the  cost  of  a slight  increase 
in  F- independent  work.  For  the  record,  here  is  such  a subprogram.  Note  that  the  input  in- 
formation in  F is  to  be  arranged  differently,  too. 

SUBROUTINE  POLSCN  (X,  F,  N) 

C INPUT  MUST  SATISFY  THE  FOLLOWING. 

C IF  X(I-l)  .NE.  X(I)  = X(I+J)  .NE.  X(I+J+1),  THEN 

C X(I+L)  = X(I)  AND  F(I+L)  = (D**  ( J-L)  ) F{X  ( I)  ) , L=0 J. 

C (HERE,  BY  DEFINITION,  X(0),  X(N+1)  .NE.  X (I ) , 1= 1 , . . . , N. ) 

DIMENSION  X(N),F(N) 

NMl  = N-1 

IF  (NMl  .LE.  0)  RETURN 

DO  3 NEXTP1=2,N 

IF  (X(NEXTPl)  .NE.  X(l))  GO  TO  4 

3 CONTINUE 
NEXTPl  = N+1 

4 DO  10  K=1,NM1 

NEXT  = NEXTPl- 1 
FLOATK  = FLOAT (K) 

NMK  = N-K 
DO  9 1=1,  NMK 

IF  (NEXT  .EQ.  I)  GO  TO  5 

F(I)  = F(I)/FLOATK 
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GO  TO  9 


NEXT  = NEXT+1 


IF  (NEXT  .GT.  NMK)  GO  TO  7 

IF  (X(NEXT+K)  .EO.  X(MEXT))  GO  TO  5 
7 F(I)  = (FINEXT)  - F(I) )/(X(I+K)  - X(I)) 

9 CONTINUE 

10  NEXTPl  = MAXO (2,NEXTP1-1) 

RETURN 


END 

We  do  not  bother  to  carry  out  here  the  extension  of  this  routine  because  it  is  straight- 
forward. Aside  from  an  initial  transfer  of  F(i,j)  to  D(j,i),  all  i.j.  only  two  state- 
ments , 

F(I)  = F(I)/FLOATK 

F(I)  = (F(NEXT)  - Fd)  )/(X(I+K)  - X(I)) 


need  to  be  put  into  a loop  over  the  data  sets,  with  the  difference  X(l+K)  - X{I)  formed 


outside  such  a loop  (and,  of  course,  F replaced  by  D(j,.)  ) . 

We  close  with  an  example  in  which  the  "matrix"  B is  three-dimensional,  taking  vectors 
to  matrices,  viz.  complete  cubic  spline  interpolation.  A typical  implementation  of  this 
scheme  (see,  e.g.,  [2])  starts  off  with  an  array  C,  of  dimension  (4,n+l) , which  contains 
the  following  information  initially: 

C(l,i)  = f (x^) , i=l,...,n+l 


C(2,l)  = f'(x,),  C(2,n+1)  = f'(x  . 

i n+ 1 

This  says  that  the  data  (X^f)  about  f in  this  scheme  consist  of  the  vector 

( f (x, ),..., f (x  ,),  f'(x,),  f'(x  ,)).  After  passing  through  a subroutine 
1 n+1  1 n+1 


SPLINE  (X,  C,  N) 

the  array  C contains  the  coefficients  of  the  polynomial  pieces  which  make  up  the  inter- 
polating cubic  spline,  i.e., 

C(j,i)  = g*^  (Xj)/( j-1)  ! , j=l,...,4  and  i=l,...,n 

For  an  extended  version,  it  would  seem  reasonable  to  introduce  a separate  input  array, 
F say,  with 
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(F(l) F(n+3))  = (f(x,) f(x  f'(x,),  f'(x  ) 

1 n+1  1 n+1 

The  calling  statement  of  the  extended  version  then  might  be 

SPLNEE  (X,  F,  N+3,  M,  C,  N) 

with  F and  C dimensioned  internally  as 


F(N+3,M),  C{M,4,N) 


Thus,  if  SPLNEE  is  used  as 

Consequently,  bicubic  spline 
out  by 


SUB  I in  the  corollary  above, 

n.  = N+3  , r.  = 4*N 

1 1 

interpolation,  on  a mesh  (x^) 


then 


n+1 

1 


by 


(y^) 


m+l 

1 


would  be  carried 


CALL 

SPLNEE 

(X, 

F, 

n+3 , 

m+3 , 

c. 

n) 

CALL 

SPLNEE 

(V, 

c. 

in+3 , 

4*n, 

F, 

m) 

with  F initially  of  dimension  (n+3,  m+3)  and  containing  the  data 


f(x^,y,) 

■■  f'^l-^m+l’ 

f(x^+l,yi)  .. 

••  ^'%+l'W> 

%'’'n+r^ 

’ ^'^n+l'^m+l’ 

fx(xi,yi) 

^x*''l'^m+l) 

f (x  ,y  ) 
xy  11 

f (*1 »y  ) 

xy  1 m+l 

^x^^n+l'^l’  •• 

••  ^x^^+l'^m+l 

)f  (X  ,y, 

xy  n+1  1 

)f  . 1 'Y  +1  * 

xy  n+1  m+l 

After  the  two  calls,  F contains  the  polynomial  coefficients  of  the  interpolating 
bicubic  spline, 

F(i+l,r,  j+l,s)  = 0/8x)^0/3y)^g(x^,y^)  , i,j=0 3 

(7) 

r=l , . . . ,n,  s=l , . . . ,m 

Note  the  difference  between  this  way  of  storing  the  coefficients  and  the  customary  way 
followed  by  the  various  available  routines  which  return  the  coefficients  in  some  array  COEP 
containing 
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COEF(i, j,r,s)  = F(i,r,j,s) 


The  coefficient  array  F,  organized  as  in  (7),  lends  itself  easily  to  evaluation  by  <xtf;id- 
ed  univariate  evaluation  routines. 

In  summary,  the  approach  to  tensor  products  advocated  here  allows  one  to  do  the  detail- 
ed programming  work  in  the  univariate  context.  The  resulting  programs  are  then  strung 
together  to  give  or  evaluate  a tensor  product  interpolant  (or,  effect  multiplication  by  a 
tensor  or  Kronecker  product  of  matrices)  with  an  ease  which  mirrors  the  ease  of  the  mathema- 
tical construction  of  tensor  products. 
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